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Abstract 

A new version of the Graeffe algorithm for finding all the roots 
of univariate complex polynomials is proposed. It is obtained from 
the classical algorithm by a process analogous to renormalization of 
dynamical systems. 

This iteration is called Renormalized Graeffe Iteration. It is glob- 
ally convergent, with probability 1. All quantities involved in the 
computation are bounded, once the initial polynomial is given (with 
probability 1). This implies remarkable stability properties for the new 
algorithm, thus overcoming known limitations of the classical Graeffe 
algorithm. 

If we start with a degree-d polynomial, each renormalized Graeffe 
iteration costs 0(d 2 ) arithmetic operations, with memory 0(d). 

A probabilistic global complexity bound is given. The case of uni- 
variate real polynomials is briefly discussed. 

A numerical implementation of the algorithm presented herein al- 
lowed us to solve random polynomials of degree up to 1000. 
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1 Introduction 



Graeffe Iteration was one of the most prestigious XIX century algorithm for 
finding roots of polynomials. At that time, computations were performed 
by hand, by people payed specifically to perform those computations. They 
were called calculateurs or computers [ [21H . 

Let / be a univariate polynomial, of degree d. Its Graeffe iterate is defined 

as: 

Gf(x) = (-l) d /(v^)/(-v^) (1) 

This defines a many-to-one mapping in the space of all degree-c? polyno- 
mials (real or complex, as wish) . The effect of this mapping is to square each 
root of /. 

After a few Graeffe iterations, the roots of G k f have (hopefully) incom- 
mensurate moduli. This is not true for complex-conjugate roots, which can 
be worked out in a different way. 

Assume, for simplicity, that / is a complex polynomial, with no two roots 
of the same modulus. Then, G k f can be written 

d 

= £a?V (2) 

i=0 

where af^ is given by the (d — i)-th symmetric function of the roots of G k f . 
Therefore, a\ is dominated by: 

{-l) d -Xi 2 \2 2k ■■■Cd-i* (3) 

where 0, ... , Q are the roots of / ordered with decreasing modulus. (A 
more rigorous statement of this will appear in Section [5]) 

Therefore, —a^/a^ is a good approximation for ( d _ i+ i 2 . 

Hence it is computationally easy to approximate \Q\ for all i. Although we 
also obtain argQ mod 2 1 ~ k ir, we will discard this information in this paper 
to avoid additional complications. There are many classical algorithms to 
recover the actual value of Q. See Pan for a discussion. [] 

In this note, we apply a suitable non-uniform change of coordinates 
(renormalization) to the Graeffe iteration operator, to make it 'convergent' 
with probability 1. 

1 Added in the revised version: We deal with this issue in [B8J. 
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The algorithm obtained by this change of coordinates will be called Renor- 
malized Graeffe Iteration. We use the following systems of coordinates for 
each iterate G k f of the Graeffe method applied to /: The coefficients of G k (f) 
and all related intermediate computations will be represented in scaled po- 
lar coordinates, where a complex number w is represented by 'magnitude' 
2~ fc log 2 \w\ and 'argument' argw 6 [— 11,71]. Calculations will always be per- 
formed 'in coordinates'. The 'magnitude' variables of G k f 'in coordinates' 
will converge with probability 1 (Theorem [1] below). 

The precise construction of the Renormalized Graeffe Operator is post- 
poned to Section ^. 

We also claim that the Renormalized Graeffe algorithm compares well 
with available numerical software or theoretical algorithms. 

However, in this paper, we are considering a modified problem: our al- 
gorithm is designed to find the absolute values of the roots, not the actual 
roots. Therefore we will compare its complexity to the complexity of finding 
the absolute value of the roots by other existing algorithms. 

In [2!|, we explain how to modify this algorithm to obtain the actual 



roots, without endangering the complexity estimates. 

The algorithm presented here has arithmetic complexity O(ol 2 ) for each 
iteration, and memory size 0(d), where d is the degree of the polynomial. 

The number of iterations will be bounded also in terms of a probability of 
failure (Theorem [| below). The authors personally believe that this is only 
possible due to the clean mathematical structure of Graeffe iteration. Our 
bound improves previous probabilistic bounds (in the sense of probability of 
success) on the complexity of solving polynomials. (See Renegar and 
S hub- S male |f45| ). 



Also, our algorithm compares well with practical software, like for in- 
stance the algorithm in Matlab (running time 0(d 3 ) and memory 0(d 2 )). 
Instead, Renormalized Graeffe Iteration seems to run in time 0(d 2 ). It also 
seems much more stable, for most (in a probabilistic sense) complex poly- 
nomials of degree, say, 1000. (See Section for a discussion of preliminary 
experimental results). 

In order to compare with deterministic algorithms one should bear in 
mind that the algorithm presented here is probabilistic and may be quite 
slow on a set of non-zero measure. Also, the complexity of deterministic 
algorithms such as |3^, f|8, 49j can be given in terms of number of arith- 



metic operations or in number of bit operations. We believe that estimates 
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of the order of 0(d 1+a ), a > in the number of operations are made at the 
expense of a large increase in the working precision, of the order of 0(d 2+a ) 
bits. However, this is not a lower bound and it may as well happen that 
those algorithms can work with substantially smaller precision on a large set 
of inputs. 

In comparison, our Renormalized Graeffe Iteration was designed for floating- 
point arithmetic (our results assume an arbitrary, but fixed floating-point 
precision). However, the details of the rounding-off analysis of renormalized 
Graeffe iteration will not be discussed in this paper. 

Instead, some preliminary numerical results are presented in Section |6|. 
They support the empirical fact that typical random degree 1000 polynomi- 
als can be solved within a precision of 64 bits of mantissa (IEEE 854 long 
double). We expect a factor of a polynomial in log<i bits of mantissa to be 
sufficient in general, with probability 1. (This is a conjecture). 

The result in Theorem [I] holds for "reasonable" probability distributions 
on the space of polynomials. By reasonable, we mean all probabilities with 
bounded Radon-Nikodym derivative with respect to Lebesgue probability in 
the projectivization of coefficient space. 

Theorem 1. There is a renormalization of the Graeffe iteration, such that 
if f is a degree-d polynomial (in a measure theoretical sense) then with 
probability 1 this renormalized Graeffe iteration produces d + 1 sequences, 
each one converging to some hi, s.t. log \Q\ = hi — h i+ i, and d, ■ ■ ■ , Q are 
roots of f . Moreover, each iteration can be performed in 0(d 2 ) arithmetical 
operations and all iterations can be performed with memory 0(d). 

This theorem is constructive, in the sense that an explicit construction of 
the renormalized Graeffe iteration will be given. 

In Section |2|, we discuss the precise meaning of renormalization in our 
context. Its main consequence, will be to produce an algorithm operating on 
a bounded set of numbers. This solves the main stability problem of classical 
Graeffe iteration that prevented it from finding all roots at once. See for 
example Henrici's comments on FFT-based Graeffe iteration |18| (Vol III, 
last paragraph of p. 69). The definition of renormalization will outlaw 
FFT-based Graeffe iteration. Indeed, although FFT is known to be stable 
with respect to vector norms, it is not component-wise stable with respect 
to the relative error. This means that some of the coefficients of the k-th 
Graeffe iterate may have a large relative error, and hence some of the roots 
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will be extremely inaccurate. This may be disastrous if one wants to retrieve 
all the roots at the same time. 



Theorem 2. Let f be a random complex polynomial of degree d. Let b > 
1 + log 2 d. Then, with probability 1 — 5, k steps of the Renormalized Graeffe 
Iteration will approximate the log \ Ci\'s with relative precision 2~ b , where 

k>c 1 + c 2 log 2 b - c 3 log 2 5 

and c-i, C3 are universal constants. The constant c\ depends on the choice of 
the probability distribution, and on d. 

Whenever speaking of random polynomials, we like to consider the nor- 
mally invariant probability density introduced by Kostlan (See also Sec- 
tion f|). However, the above mentioned result is true for any reasonable 
probability distribution. 

The experimental results in figure ^| support the conjecture that, under 
Kostlan's probability distribution, we can fix 

ci = c 4 log 2 d 

where C4 ~ 2 is a universal constant. 

We will briefly discuss the real case, and how to deal with complex- 
conjugate roots or roots with same modulus in Section |5|. 

Historical remarks. The Graeffe iteration was developed independently 
by Dandelin (1826), by Graeffe (1837) and by Lobachevsky (1834). We call it 
Graeffe iteration to conform with most of the literature. See Householder |2(| 



for early references and priority questions. See Dedieu || for an application 
of Graeffe 's algorithm. 

Important theoretical results were obtained by Ostrowskii [33J] in 1940. 
Also, by that time, numerical analysis books mentioned Graeffe iteration as 
the preferred algorithm for zero-finding (see e.g. Uspensky [^TJ Page 318. 
For another early computer implemented algorithm, see Bareiss |l|, 0] and 
also Blish and Curry ||). 

With the advent of digital computing, the practical use of Graeffe itera- 
tion seems to have been forgotten. 

Most popular zero-finding algorithms seem to be based now on QR iter- 
ation (Matlab) or in a several steps, root-finding plus deflation scheme, (e.g. 
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Jenkins and Traub 



See, however, Cardinal 



Edelman and Mu- 



rakami [13|], Emiris, Galligo and Lombardi ||14|| , and Toh and Trefethen |50 

In a more theoretical perspective, Graeffe iteration is considered as a sort 
of pre-conditioning for polynomial splitting. Splitting a polynomial means 
factorizing it into one factor with large roots, and another with small 
roots. Splitting is used to obtain extremely fast theoretical algorithms (see 
Schdnhage |j| gj|, Kirrinnis |3|, Neff and Reif H, Bini and Pan §, Mour- 
rain and Pan 0, Pan f3|, 0, g|, Pan et al. |9f 



Bini and Pan ^ 
Pan |36|, 0, Pan et al. p9| , Malajovich and Zubelli 
The main practical difficulty for those algorithms seems to be the large 
precision required by Graeffe iteration. Also, those algorithms are quite close 
to known lower bounds on topological complexity (See Vassiliev |52|). For 
a related lower bound see Novak and Wozniakowski 



An important paper by Grau in 1963 |17j laid some of the bases for a 
version of Graeffe iteration adapted to digital computers. He identified the 
problem of the increasing numerical range. During Graeffe iteration, some 
of the coefficients can become so large that the floating point system cannot 
accommodate them anymore. 

While most of the literature suggests to find one root at a time and then 
use deflation, disregarding some stability problems (See e.g. Henrici [0]), 
Grau proposed a globally convergent algorithm. Grau's algorithm would 
involve only bounded quantities. 

As far as we know, that paper was completely forgotten. The algorithm 
suggested by Grau has complexity 0(d 2 ) and memory usage of 0(d 2 ). It 
may be considered as the precursor of the one we shall introduce below. 



2 Iterative Algorithms and Renormalization 



In this paper, we will produce a version of Graeffe iteration that has bounded 
numerical range, for most input polynomials. The crucial concept in the 
construction of this algorithm is the idea of renormalization. 

Renormalization is a tool used in understanding the qualitative behav- 
ior of iterative phenomena that range over different scales. A rich the- 
ory of renormalization exists for one-dimensional dynamical systems. See 
Feigenbaum McMullen [2Jj, and De Melo-Strien 

HI- 



As for the multi- 



dimensional case see Palis- Takens 

We shall illustrate what we mean by renormalization in the setting of 
iterative algorithms by a well-known example and then proceed with the 
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definitions 



Example 1. Let / : C — > C be a smooth function. Newton iteration 
associates to a point x the point N(x) = x — f(x)/f'(x). The sequence of 
Newton iterates converges to a zero of /, provided that x is picked close 
enough to a non-degenerate zero. 

We assume now that f(x) = + f'(0)x + h.o.t. , and that f(0) ^ 0. We 
will consider now the renormalization of the iterative algorithm (or mapping) 
x i— > N(x). Although we are dealing here with a simple example (the answer 
is always 0), its renormalization will have some of the main features of the 
renormalized Graeffe iteration, yet to be defined. 

The basic idea, therefore, is to look at Newton iteration with a variable 
microscope of variable lens. We look at the mapping: 

N e = h e -2 oNoh e (4) 

where h e means the homothety of ratio e. 

When e tends to zero, N t tends to the map: y h- > r yy 2 , defined from the 
disk D(|7 _1 |) of radius |7 _1 | onto itself. Here, 7 is a parameter chosen equal 
to / (2) (0)/2/'(0). (Proof of this fact: Taylor's Theorem. The choice of the 
radius I7I -1 makes the limiting map surjective). 

This process (which we call renormalization) gives us qualitative infor- 
mation on the dynamics of Newton iteration near a root. We can summarize 
this information in an eventually commutative diagram (commutative in the 
limit): 



£>(|7 _1 |e) 3 x 



N 



£>(| 7 



-H e -2; 



3N{x) 



y £D(\r 



(5) 



72/ 2 e£>(|7 



-II 



In the example above, the homothety {h e ) _1 is generating the renormal- 
ization group. In general, we want to consider renormalizations that lead to 
a commutative diagram, on the limit: 

R 



y 



G 



G{y) 



R? 



R{y) 

G R 

R 2 (G(y)) 



(6) 
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Above, G is the original algorithm (possibly after some uniform change 
of coordinates). We denote by R the renormalization map, and by G R the 
renormalized version of our iteration. Usually, G R depends on a renormaliza- 
tion parameter. However, we want G R to converge to a limiting map, with 
a simple dynamics. Moreover, computation of G R should be 'stable', in a 
precise sense. 

This suggests the following heuristics, in the case of Graeffe iteration: 
we would like to consider a polynomial / represented by its roots. More 
precisely, a polynomial / can be represented by the vector (log arg(£j)) 
eK^x T d , where Q is the i-th root of /, and T denotes the additive group 
]R/27rZ. In that case, Graeffe iteration is just multiplication by 2. 

Renormalization would be a division by 2 of the log of the radii of the 
roots. We would have: 

l d xft (log|Ci|,argCi) — (r,9) eR d x T d 

(r,0)H-»(r,20) 

R d x T d 9(2 log |0|,2 argO mod 2vr) ► (r, 26 mod 2n)eR d x T d 

R k+i 

(7) 

Therefore, the limit of the renormalized map should be the map: (r, 6) \— > 
(r,28). Of course, we do not know in advance the roots of the polynomial. 
Instead, we will produce a chain of commutative diagrams 'converging' to 
diagram (^). In order to do that, let us assume that a polynomial 

a + a\X + a 2 x 2 + ■ ■ ■ + ddX d (8) 

is represented by the vector 

(ai, at) = (log |di| , arg(ai)) (9) 

In this paper, for clarity of exposition, we will ignore the case where some 
of the a, is zero. We can do that because most of our results hold 'with 
probability 1'. In practice, polynomials with a zero coefficient do arise. See 
Section ^ for implementation comments. 

If the operator G represents Graeffe iteration in this new system of coor- 



x2 
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dinates, we will have: 



Dd+l x Jd+1 3 Qk(f) 

G 

i d+1 x T d+1 3G k+1 {f) 



R k 



* (a k , a k ) eR d+1 x T d+l 



Rk +i 



(Ok) (10) 
- (a k+ \a k+1 )eR d+1 x T d+1 



Above, a k = (Sq, • • • , d^) and a k = (a k , • • • , a k ). Also, A; is a superscript, 
not an exponent. However, R k denotes the k-th iterate of R. 

In this diagram, R maps (r,9) into (r/2,9). Although this defines Gk 
as a mapping, the algorithmic construction of Gk is postponed to Section 
[| On the limit, Gk 'converges' to the mapping: (a, 9) i— > (a, 20). This 
will ensure that lima fc exists. We will show that this limit satisfies with 
probability 1 

lim d k - lim a k +l = log \Cd-j\ 

fc^oo fc^oo 

where Cd-j is the (d — j)-th root of the original polynomial, the roots being 
ordered by decreasing modulus.. 

Also, by using the renormalized Graeffe iteration Gk instead of the clas- 
sical Graeffe iteration G, we will be able to bound all the intermediate cal- 
culations to a compact, depending on the input /. This will imply several 
stability results. 

In order to define formally what we mean by a renormalized algorithm, we 
need to state the conditions that we expect Gk to satisfy. These conditions 
will define a class of algorithms, that we call renormalized iterative algo- 
rithms. Some examples: usual Newton, renormalized Newton, all reasonable 
algorithms based on a contraction principle and, of course, Renormalized 
Graeffe Iteration (To be constructed). 

Before defining what we mean by a renormalized algorithm, we will briefly 
discuss the notion of algorithm. There are many possible definitions (See 
Blum, Cucker, Shub and Smale [0]). 



10 



Definition 1. An iterative algorithm M is a Blum-Shub-Smale (BSS) ma- 
chine over R, modified as follows: 

1. If the input is in R l x T m for a pair (Z, m) , then the output is in R l x T™. 
The integer / + m is called the input size. If the input is denoted 
by x, the output is denoted M(x), and M k means the composition 
M o M o • • • oM 

Also, we say that the iterative algorithm 'computes' the function x \— > 
lim^oo M fc (x) 

2. Computation nodes are allowed to perform the following elementary 
functions: polynomial evaluation, absolute value, (real) logarithms, 
(real) exponential, sine, cosine, rational functions, and compositions 
of those. 

3. M e (x) will denote the result of approximating M(x) by allowing each 
operation with non-integer parameters to be performed with relative 
precision e. (This is also known in numerical analysis as the (1 + e) 
property.) The parameter e is allowed to vary in (0, |). The machine 
is supposed to "know" e. This means that the value of e can be used 
in intermediate calculations. Also, the approximation is performed by 
some prescribed algorithm, (e.g., IEEE arithmetic with — log 2 e bits of 
mantissa). 

4. For all e, the arithmetic complexity of M applied to the input x is the 
number of elementary function evaluations (from the list of item 1) 
and branchings performed with input x. We require the arithmetic 
complexity of the approximate algorithm M e applied to the input x 
to be the same as the arithmetic complexity of the original algorithm 
M applied to input x. 

Remark 1. Item 1 will allow us to distinguish between real and angular 
variables, the latter one being defined modulo 2ir. This will simplify notation 
when we speak of the distance between points in M 1 x T m . Let d be that 
distance. 

Remark 2. It has been argued that the outcome of a branching node would 
not be well-defined in the presence of numerical error. This is not true in the 
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definition above. The branching nodes of the machine M can be assumed, 
without loss of generality, to branch on queries of the form y > or y > 
or y = 0. When the machine M gets replaced by M e , the branching nodes 
stay formally the same. The value of y, however, is contaminated with a 
certain numerical error. It is still a perfectly defined real number, and it may 
be compared to zero. Thus, the branching nodes branch 'correctly', for a 
slightly perturbed input. 

This would lead to disastrous results if an approximate machine enters 
a 'loop' due to numerical errors. Item 4 in the definition above is there to 
preclude that sort of loop, by ensuring that the arithmetic complexity of each 
iteration does not depend on the working precision. 

Definition 2. A function tp can be computed in finite time if and only if 
there is a BSS machine over R modified as above items 2, 3 and 4, that 
computes ip. 

One consequence of the previous definition is the following: suppose a 
certain function tp can be computed in finite time. Then, its branching set 
is given by a finite set of equations. Outside the branching set, ip can be 
written (locally) as a composition of elementary functions. Moreover, only a 
finite number of such compositions may appear, one corresponding to each 
set of possible branchings. 

A few definitions are in order now. In the sequel we will need to use 
algorithms depending effectively on a parameter k G N. Those will be 
given by a machine M with two inputs, say k and /. However, we will 
denote the output as M k (f) and we will write M k for each of the input- 
output mappings obtained by restricting that machine to some fixed value of 
k. Also, in that situation we will speak explicitly of the algorithm (M k ) keN 
or M k for short. 

The sequence (Mfc)^ =1 can be considered as a sequence of mappings from 
M. 1 x T m into itself. The orbit of / by the sequence (M fe )^ =1 is the set 

orb/ ^ (/, Mi(/), M 2 o Mx(f), . . . ) C M ; x F . 

This set should be understood as the orbit of / by the non-autonomous 
dynamical system M k (Not the semi-group !) The closure of orb/ will be 
denoted by orb/. 
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We shall say that a subset of M 1 x T m has full-measure if its complement 
is contained in a set of null measure. We say that some property is true 
almost everywhere (a.e.) if that property is true in a full-measure set. We 
can now define a renormalized algorithm and make precise our concept of 
renormalization: 

Definition 3. The algorithm (M k ) ke ^ (or M k for short) is said to be a renor- 
malized iterative algorithm to compute ip : M} x T m — > W, f \— > tp(f), where 
< s < I, if, and only if, it satisfies the Axioms 1 through 4 below. 

Axiom 1 (Consistency). For almost every / G M 1 x T m we have 
lim (tt o M k o M fc _! o ■ ■ ■ o M 1 ) (/) = 

fc— +oo 

where 7r is the projection of R l x T m onto the first s coordinates of R l . 

Axiom 2 (Arithmetic Complexity). The arithmetic complexity of M k 
with input / is bounded in terms of the size I + m of /, independently of k 
and the coefficients of /. 

Axiom 3 (Propagation). For almost every / G 1R ; x T m , there exists a 
compact neighborhood V C M 1 x T m of orb/ (under {M^}) and C such that 
M k \ v is Lipschitz with constant C k and eventually C k < C . 

Axiom 4 (Stability). For almost every /eM'x T m , there exists a com- 
pact neighborhood V C M z x T™ of orb/ (under {M k }) and 5 G R such that 
Wg <E V and Vfc, 

d(M kte (g), M k {g))<eB . 

Our concept of renormalized iterative algorithm subsumes several 'rea- 
sonable' properties of iterative algorithms. Axiom 1 allows our algorithm to 
carry more information than what is actually required at output. Yet, we 
want our algorithm to produce a sequence converging to the expected result, 
for almost every input. Axioms 3 and 4 will rule out unstable algorithms. 
The idea behind Axiom 2 is that any honest iterative algorithm should have 
bounded arithmetic complexity for each iteration. This prevents the use of 
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multiple precision arithmetic to obtain stability at the expense of (possibly 
exponentially many) extra arithmetic operations. 

We shall now explore some consequences of the definition of renormalized 
iterative algorithm we proposed above. More specifically, our first goal is 
to obtain a (non-uniform) error bound on the result of iterating k times a 
renormalized algorithm with precision e. 

Lemma 1. Let M be a renormalized iterative algorithm . Then, for almost 
every f and for each k, we have that for sufficiently small e > 

d (M k , e o.-.o M lje (/), M k o ■ ■ ■ o M 1 (f)) < kA k Bt . 

Here, A and B, depend on f , but not on k. In particular, that will be true 
for values of e of the form 

< P 
e ~ kA k B ' 

where p is defined below. 

As an immediate consequence of Lemma [l] we have that 

|| (io%o...o M 1>e )(f) - (vr o M k o ■ ■ ■ o Mi)(/)) || a < kA k Be . 

Proof of Lemma [J. We now describe the construction of the constants in the 
statement of Lemma |I[ Combining Axioms 3 and 4, there exists a compact 
neighborhood W of orb/ such that: 

1. Every mapping M k is Lipschitz on W. Moreover, there exists a constant 
A such that for every k, the norm of the Lipschitz constant of M k is 
uniformly bounded by A = max (1, sup fc C k ), C k as in Axiom 3. 

2. There exists 6q such that Ve < eo we have 

d(M k>t (g),M k (g))<eB , Vg G W . 

We then define p as the minimum of eo and the distance of orb/ to the 
boundary of W. 
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We may now conclude the proof of Lemma [I]. Set 

</ = (M fc _ 1 o...oM 1 )(/) , 

and 

/i=(M fc _ 1>e o...oM 1>e )(/) . 
In that case, we want to bound 

d(M k , e (h), M k (g)) < d(M k>e (h), M k (h)) + d(M k (h), M k (g)) 
If h G W, one can bound 

d(M k (h),M kje (h))<Be 



and 



By induction, 



and hence 



The condition 



d(M k {h),M k {g))<Ad(g,h) . 

d(g,h) < (k-l)A k - l Bt , 
d(M k , e (h),M k (g))<kA k Be . 



e< P 



kA k B 

guarantees that h G W. □ 

If we use a Turing machine model (or any other classical discrete com- 
plexity model), we can perform all the operations of M in finite precision, 
and obtain: 

Lemma 2. Let M be a renormalized iterative algorithm. For almost every 
f , assume that the truncation error is bounded by 

|| (7T0M fc 0...0M 1 )(/)-^(/)|| 2 <^ 2fc , 

where E = E(f) G (0, 1). Then, the complexity of approximating <f(f) with 
precision S is 0((log 2 |) 1+Q ) where a > is arbitrarily small. 

15 



It is assumed above that the cost of arithmetic with I bits of mantissa is 
0(l 1+a ), for all a > 0. See J3J pp 78-79 for a sharper bound on the complexity 
of long integer multiplication. 



Proof of Lemma Let 



so that 



log 2 



log 2 E 



E 2k <-. 



Choose e so that Lemma || holds, i.e., e < p/kA k B, and such that 



kA Be < 5/2. This can be done for 



e < Ci<5(log<5 



for some constants ci and c 2 dependent of /. The total cost of computing 
is O (a(log 2 e _1 ) 1+a ), where a is the arithmetic complexity of the algorithm 
and a is arbitrarily small. Since / is fixed, a is a constant and k < 
log 2 log 2 (5~ 1 the total cost is 



^((log.r 1 ) 1 ^) 



for all a arbitrarily small. 



□ 



The complexity bound of Lemma [2] is non-uniform in /. It is also non- 
effective, in the sense that we give no procedure to estimate k and e without 
the knowledge of p, B and A. Indeed, those quantities may depend on /. It 
would be of some interest to bound those quantities in a probabilistic setting, 
similar to Theorem p|. 

Definition 4. Let M be an iterative algorithm to compute $ : F i— > <&(F), 
i.e., lim^oo M n (F) = $(F)). We say that is a renormalization of M if 

1. Mfc is a renormalized iterative algorithm to compute f(f)- 
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2. There are functions ip and 77, defined almost everywhere and com- 
putable in finite time, such that the diagram 



^— (p(ip(F)) 

commutes. 

3. There is a function R, computable in finite time, so that the diagrams: 



M 



R k (^(F)) 



M(F) Rk+1 °% R k+1 (ip(M(F))) 



commute. 



Theorem [I] can be stated now in a more concise way: Let ( : / £(/) 
be the function that associates, to any univariate degree d polynomial /, its 
roots Cij • • • j Cd ordered by decreasing modulus. We have the following: 

Theorem. There is a renormalization of the Graeffe iteration to com- 
pute |C(/)| = (|Cl(/)l> • • • > 1Cd(/)D- Each iteration has arithmetic complexity 
0(d 2 ) and needs memory 0(d). 



3 Recurrence Relations and the Renormal- 
ization of Graeffe 

It is time to construct the renormalized Graeffe iteration. Let f(x) = /o + 
f\x + ■ • • + fdX d . Let h = Gf be its Graeffe iterate. The coefficients of h can 
be written as : 

0<i-j<d 
0<i+j<d 
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For convenience, we rewrite QTTj ) as: 

hi = (-l) d+ 7, 2 + 2 £ i -IV'" 7, .;/,-.; (12) 

1<j <min(2,<i — i) 

The next step is to write those equations in terms of the log of the coeffi- 
cients. More precisely, we will have to deal with the following two quantities: 

/i° g = log|/,| (13) 



r* = arg/i (14) 



It is possible now to construct the renormalized Graeffe iteration Gj.. 
Recall from diagram (pTDj) that this iteration will map 2~ k f log and / arg into 



2- k - 1 h log and /i arg 

For that purpose, we introduce the notation: 



jk def (2~^ y'°s / arg ) 



h k+1 = (2" fc - 1 /i log ,/i arg ) 
We also introduce operators: 



(15) 



(x,a)®(y,fl) = (x + y,a + P) 
(x, a) = f (Ax, Aa) 



and 



(x,a)ffl fc (y,/3) 



def 



fx + 2 'log b|, a + argz) 



(16) 
(17) 
(18) 



e m+2 fc a' + e i/3+2 fc j/ 



. u;ft ( e to+2*» + e ^+2*v^ ( 19 ) 



We remark that the purpose of the sub-index k in the above formulae 
is to keep track of the degree of the renormalization. For operations which 
do not change, we omitted the sub-index. The operator \tl stands for the 
multiplication of a renormalized value by a (non-renormalized) constant z 
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Also, binary renormalized operations are defined for operands with the 
same renormalization index. Therefore, one should first convert f k to f k+1 
before attempting to 'multiply' it with a factor of renormalization index of 
order k + 1. This conversion will be implicit in the formulae below. 

Equation flT2] ) becomes: 




(20) 

where = (— Recall that above, k is a superscript, not an expo- 
nent. The 'renormalized operations' above are easy to implement in terms 
of the classical ones. The most delicate being the renormalized sum, so we 
give here our preferred algorithm 

Example 2. How to compute the 'Renormalized sum': 

(c, 7 ) := (a,a)m k (b,P) 



If a > b, do : 

s = exp ia + exp (i(3 + 2 k {b — a)) 
c = a + 2~ k In \s\ 
7 = arg(s) 

else 

s = exp i/3 + exp (ia + 2 k (a — b)) 
c = b + 2~ k In \s\ 
7 = arg(s) 

endif 



We remark that in the above formula, the complex arithmetic operations 
can be performed in terms of real elementary ones. Moreover, in numerical 
implementations if k is large enough, as compared to e, it may be faster to 
approximate (0,7) with (a, a) or (b,/3), whichever is larger. 

In order to finish the proof of Theorem 0, we still need a few remarks about 
the renormalized Graeffe iteration, which we just constructed. Clearly, in 
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order to compute formula (|2"0|), we only need memory space of the order 0(d) 
and time of the order 0(d 2 ). 

When k — > oo, the quantities |/*| are all convergent. 

If Gk is the renormalized Graeffe iteration, we define Goo as the limit 
of Gk when k goes to infinity. This means that renormalized sums EE^ are 
replaced by their limit EBoo, where: 

{(a, a) if a > b 

(b,(3) ifa<6 
undefined if a = b 

It is easy to see that Gk — > G^ almost everywhere, pointwise and in 
the C l topology. In order to avoid a rather tedious calculation, we can 
establish this fact from the pointwise C 1 convergence almost everywhere 
of (a, a)S k (b, (5) to (a, a)ffloo(6, f3) and similarly for the other renormalized 
operations. 

We define: 

tp(fo, ■ ■ ■ , fd) = (log |/o|, - - - ,log|/d|; arg |/o|, - - - , arg|/ d |) 
r}(r , • ■ ■ ,r d ;a , • ■ ■ ,a d ) = (exp (r - n), ■ ■ ■ , exp (r d _i - r d )) 

R(r , ■ ■ ■ , r d ] oo, • • • , ad) = ( -pr, • • • , 7?; ao, • • ■ , a>d) 



We define $ as the function that associates to a polynomial / the values 
|Ci|, • ■ ■ , |Cd| where Q are the roots of /, ordered by decreasing modulus. 
Then, we set 

cp = r]^ 1 o $ o ip -1 



With the definitions above, Gk is indeed a renormalization of the classical 
Graeffe algorithm to compute $: 

Proposition 1. Renormalized Graeffe Iteration is a renormalized iterative 
algorithm (in the sense of Definition^) to compute $. 
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The proof of this proposition will conclude the proof of Theorem |T|. In 
order to establish Proposition [I], we should check that our algorithm, as 
described in equation (^), satisfies Axioms 1 to 4. 

Axiom 1 is verified by construction. Axiom 2 follows from the recurrence 
formula ©. 

The proof that our algorithm satisfies Axiom 3 will require a technical 
lemma. Before stating it, we recall some notation: 

orb(/) = {/; G a (/); G 2 o G x (f)- . . . ; G k o - • • o G 1 ( f);...} is the orbit of 
/ under the sequence (Gk)- Its closure is denoted by orb(/). 

We will show the following lemma, which implies satisfaction of Axiom 3. 

Lemma 3. For almost every f and any 5 > 0, there exist a compact neigh- 
borhood W C K x T m of orb(f) and an integer ko such that, Gk is a local 
diffeomorphism W — > G k (W), and such that for ko < k < oo, the derivative 
of Gk is bounded by 2 + 5 

Proof of Lemma 0. The proof is divided in several steps. 

Step 1: For any j G N, there is a full-measure set Uj such that Gj o ■ ■ • o G\ 
is well-defined, and a local diffeomorphism in Uj. Hence, there is a 
full- measure set Uj t k such that Gk ° (Gj o ■ • • o Gi) is defined on Uj t k 
and for all / G Uj t k there is a neighborhood Vf of Gj o • • • o G±(f) such 
that Gk is a diffeomorphism Vf — > Gk(Vf). 

Step 2: Let Uoo be the set of all / such that Gk is a diffeomorphism near 
(7r _1 o(^)(/) for all values of k that are large enough. Then contains 
the set of complex polynomials without roots of the same modulus. 
Hence has full measure. 

Step 3: Let U = Uoo H (f)j Uj^j fl \ f)jkUjkj- Then U has full measure. More- 
over, let / G U. Then is a local diffeomorphism with derivative of 
norm < 2 + 5/2 in an open neighborhood Vfc(g) of every g G orb/, for 
all > k (g). Indeed, if we write Gk = Goo + (Gk ~ G^), we can make 
the C 1 norm of the second term arbitrarily small, namely less than 5/2. 
We know that the norm of the derivative of Goo is precisely 2, hence 
the bound 2 + 5/2 In the particular case 5 = oo we can set ko — 1. 

Step 4: Since Gk — > pointwise in the C 1 topology and for <? almost every- 
where, we can assume that f] k> k Vk(g) contains an open ball V(g) of 
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center g, where G k is a local diffeomorphism with derivative bounded 
by 2 + 5/2. 



Step 5: Since orb(/) is compact, the union U g g orb(/) nas a finite sub-cover 
Urer v (d), and we set W = {J g& V(g). Then we set k = max 9er k (g), 
and we obtain that for any k > ko, G k is a local diffeomorphism in W, 
with derivative of norm bounded by 2 + 5. 

□ 

The proof that our algorithm satisfies Axiom 4 is divided in two parts , 
dealing (respectively) with small and large values of k. 

Lemma 4. For almost every f and for all k, there exist an open neighbor- 
hood U containing f and B > such that for all g G U and for all e small 
enough 

d(G k , e (g)-G k (g))<Be . 

Lemma 5. For almost every f , there exist k$, B > 0, and an open neigh- 
borhood U containing f , such that for all g G U and k > ko, for all e small 
enough, 

d(G k>e (g),G k (g))<Be 

Lemmas |] and |5| together imply that for almost all /, there is a neigh- 
borhood U containing / such that for all g G U, 

d(G k>e (g),G k (g))<Be 

independently of k. 

Since orb(/) is defined for almost every / and admits a compact neigh- 
borhood V, we can select a finite subcover of those neighborhoods, and hence 
find a finite B valid for all g G V. Therefore, our algorithm satisfies Axiom 
4. 

Proof of Lemma 0: Our iteration G k can be written in terms of the following 
real operations: +, — , cos, sin, arctan, multiplication by 2 k , by 2~ k , absolute 
value, exp, log. 
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The set of inputs such that a 'log of zero' or an 'absolute value of zero' 
occurs has zero measure. Therefore, for almost every /, Gk(f) is computed 
by a composition of analytic functions. Also, for almost every /, none of 
the output values is zero. 

Therefore, for every intermediate quantity xi, the derivative of any output 
y m with respect to x\ is finite (say < Am)- Therefore, a relative pertur- 
bation of e in X[ leads to a perturbation of size \xi\Di m e in y m . Thus, we set 
B = | ^2i m Di m \xi\. By continuity, a perturbation of e in each intermediate 
value xi leads to a perturbation smaller than Be, for input g in a certain 
neighborhood of /. □ 



Proof of Lemma [|: Let W be the set of all / such that G 00 (/) is well-defined. 
Recall that 

ffloo : (a, a), (b, 0) t— > (a, a) if a > b 

(6, (3) if 6 > a 

and that the operator EB^ is not defined for a = b. Therefore, W is open 
and has full measure. Let / G W and Z7 be a small connected neighbor- 
hood containing /. Let g E U, then by taking U small enough and large 
enough, we can guarantee that Gk{g) and Gk,e(g) are well-defined. Since U is 
connected, all the branching outcomes in the computation of Gk{g) and Gk, e 
are the same, hence Gk restricted to U is a composition of locally analytic 
functions. We can assume without loss of generality that all derivatives are 
bounded, hence there is a constant B such that 

d(G k>e (g),G k (g))<Be 

for e small enough, but still independent of the choice of g. 

□ 

This concludes the proof of Proposition [I] and hence of Theorem [l] 



4 Probability of Success 

Through this section, \\.\\d will denote Weyl's unitary invariant norm (|53 
III— 7, pp 137-140) in the space Vd of complex polynomials of degree at most 
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ci(See 0). If f{x) = $Xo f iX \ then 



1 



i=0 



This norm is invariant under the following action of the group U(2) of 

a (3 



unitary 2x2 matrices: if (p - 



is unitary, define 

ay + (3 s 



Thus, \\r\\ d = \\f\\. (See Ch. 12, Th. 1 of §].) 

This is in some sense the most 'natural' norm in the space of all polyno- 
mials. More information about that norm, its associated probability distri- 
bution and its applications can be found in |3|, [10], [11], |25], |26], ^7], [4^, [0 
H @- 

Let PPrf be the projectivization of normed complex vector space (Vd, \\-\\d) 
We can define the 'sine' distance in FVd by: 



dp(f,g) 



mm 

Aec 



llZ-AsH 



where £> is the usual Riemann metric in PVd- Also, there is a natural volume 
form in FVd- Normal invariant distributed polynomials in Vd correspond to 
uniformly distributed 'polynomials' in FVd- 

Let / be a random polynomial (in any of the two equivalent senses above). 
Then we may consider its roots £1, • • • ,(d as random variables. The joint 
distribution of £1, •• • , Q was studied by Kostlan in 
result below will rely on elementary estimates: 

Lemma 6. Let f be random in the sense above. The probability that 

mm — — - > 1 + e 

Kii>ioi 101 



25 . However, the 



is larger than 1 — Me, where M is a positive constant depending on d. 
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This result is also true if we chose / random with respect to any other 
probability distribution with bounded Radon-Nikodym derivative with re- 
spect to the volume form in FVd- The constant M will have to be multiplied 
by the maximum of the Radon-Nikodym derivative of the new distribution 
with respect to the volume form. 

Lemma |6| will be a consequence of the 'Condition Number Theorem' be- 
low. Let 

p(f) = min 1 — 

ICiKIOI 101 

We will interpret p(f)~ l as a condition number. Let Eg be the locus of 
ill-posed problems, i.e., the set of polynomials such that = \Q\ for some 
i ^ j. Then, 

Theorem 3 (Condition Number Theorem for Graeffe Iteration). 

dp{f, E G ) 



P(f)> 



Vd 



Therefore, the probability that p(f) > e is no less than 1 — Vol V ev qLc 
where V e ^q denotes an ev^-neighborhood. This volume is of 0{t\fd) < Me. 

Lemma 7. Let f(x) = (x — (i)g{x) G Vd, where g e V d -i- Then 



\9\\d-i 



< 



d 



i + ICi 



Proof of Lemma Q. We start with the easy case and assume that £i = 0. Set 
9(x) = Eto^*- Then, f(x) = Yl=\9i-\^ ■ Now > 

1 I |2 d I |2 (I i | |2 

II ||2 \9i\ \9i-l\ " \9i-l\ ^ M f\\7 

n»ii«-> = E = E TTrrr = E -77T * "ii/n. 

i=0 I t=l I i=l 



i J \ i — 1 / \ i 

and hence \\g\\d-i < Vd\\f\\ d . 

For the general case, we will use [7(2)-invariance of and Let 
(p be a convenient unitary matrix: 

1 r 1 Ci" 

25 



Set f v = f o ip and g v = g o tp. The choice of tp has the particularity that 
/v(0) = /(Ci) = 0. We can compute in terms of g v : 



r( y ) = VT+W y g^y) 

Using the easy case, 

||V^HCJVIU-i < VdWFWt 

By ?7(2)-invariance, 



\\g\U-i = \m\ d -i 

v 7 ^ 



Lemma 8. Let g e V d -\- Then \\g\\ d < \\g\\d-i- 
Proof of Lemma 0. 

^Wd - 2^ ~T~T\ ~ 1^ A ( A_ 1 \ s 
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Putting Lemma |7] and Lemma § together, 
Lemma 9. Lei f(x) = (x — Ci)fi'( a; ) £ ^d; where g G Vd-i- T/ien 

d -7ArF ll/lu 



□ 



□ 



26 



Proof of Theorem |3j. Let f(x) — (x — Ci)( x — ( 2 ) • • • {% — Cd) an d order the 
Cj's such that 

/<-\ 1 Ci i C2 

pit)— mm 1 - — - = 1 - — — 

iokici 101 ICil 

Define h(z) — (x — — p))(^ — C2) • • • (x — Cd) £ Eg- Then 

11/ -^11, = ||Cip(*-CO-"(*-Ci)l|d 

= [CxIpIK* - Ca) • • - (a: - Ci)IU 
< ICilP 7=L^ \\fh 



Hence, 



Hence, 



ICil 



□ 



Proof of Theorem [|. We set 5 = Me, where M is the constant such that the 
volume of an e\fd neighborhood of Eg is less than Me. With probability 
larger than 1 — 5, 

max — — ■ > 1 — e . 

ICiKIOI ICil 

We now use the fact that, for any N > (log2)/e, we have (1 — e)^ < 1/2. 
We set ki — 1 + |~log 2 e -1 ] and using N = 2 kl in the previous formula, we 
obtain: 

max Ml < (1 _ e f^ < I . 
IGKIfcl ICil 2 1 2 

An extra 1 + log 2 b iterations ensures that 

max ML < 2- 1 -" , 

ICil 2 
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provided that 

k = k\ + 1 + log 2 b . 
We claim that we have in this case that the sum 

E #-£ = C? ■■■Cf(l + e), 

i\<—<i r 

where 

|e| < 2~ b . 

Indeed, assume we have reordered the roots so that \^\ > ■ ■ ■ > It then 
follows that 

Now, for r fixed, let 
<Sk = {{h, ■ ■ ■ , i r ) \ h < ■ ■ ■ < i r and i\ + . . .i r — r(r + l)/2 = k} . 

Note that S = {(1, . . . , r)}, and S x = {(1, . . . , r - 1, r + 1)}. Also, S 2 = 
{(l,...,r-l,r + 2);(l,...,r-2,r,r + l)} 

In general, every multi-index i\ < • • • < i r may be obtained by starting 
from 1 < 2 < • • • < r and increasing one of the indices, in such a way not 
two indices are equal. Then S). is the set of multi-indices obtained after k 
steps. Therefore, we may bound #Sk+i < d #5^ to get 

J2 2-( 6+1 ) fc #S k < 2- b -\l + 2- b + 2- b d + 2~ 2b - 1 d 2 . . . ) < 2- b , 
fe>i 

under the assumption that d < 2 b . This concludes the proof of Theorem |2|. 

□ 



5 Newton Diagram Revisited 



Following Ostrowskii [p^l , we introduce Newton diagram, but now with a log 
scale. The Newton diagram of a polynomial / is the graph of the piecewise 
linear function defined by g(i) = — log|/j|. (Our construction differs from 



28 



Perfidous polynomial of degree 20, before iteration 



"per1idous.20.initial" 



Perfidous polynomial of degree 20, after iteration 



10 12 14 16 18 20 2 4 



10 12 14 16 18 20 



Generic polynomial ot degree 1000, before iteration 



50 







"genoric. 1000. initial" 
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Generic polynomial of degree 1 000, after iteratio 
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100 200 300 400 500 



700 800 900 1000 100 200 300 400 500 



Figure 1: Examples of Newton diagrams 



700 800 900 1000 



Ostrowskii's point of view. He used to define the Newton diagram as the 
convex hull of that set of points). It makes sense to renormalize Newton 
diagram, as we did with polynomials G k f. 

Therefore, we may speak of a Newton diagram on the limit. Assume (as 
we did through this point) that all roots of / have different moduli. Then 
the derivative of this diagram on the piecewise linear part over [i,i + 1] will 
correspond to log | Cd— i I ? where Q is the i-th largest root of / in modulus. It 
follows that the Newton diagram will be the graph of a convex function. 

We should consider now the general case. If / has several roots of the 
same moduli, some of the ratios will converge to the ratios of the 

coefficients of the factor of / containing those roots. The case of three roots 
of same moduli is illustrative: 

Let Ci> • • • > Cd be arranged by non-increasing moduli, and assume \0\ = 
1 0+i | = |C*+2|- Then, 



fd-i = Cl • • • Ci-1 (Ci + 0+1 + Ci+2) + • • • 
fd-i-1 — Cl • • • Ci-1 (CiCi+1 + +CiCi+2 + Ci+lCi+2) + 
fd-i-2 — Cl • • • Ci+2 + • • • 



(21) 
(22) 
(23) 
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Hence, on the limit, we approximate: 



fd-i-2 
fd-i-l 
fd-i-1 
fd—i 
fd—i 



00+1 0+2 



00+1 + +00+2 + 0+l0+2 
00+1 + +00+2 + Ci+lC»+2 
+ 0+1 + 0+2 

+ 0+1 + 0+2 



fd-i+1 

The moduli of those roots is therefore given by: 



(24) 
(25) 
(26) 



f 



d-i-2 



d-i-1 



f 



d-i-1 



fd- 



f 



d—i 



f 



d-i+1 



(27) 



The same formula extends for factors of any degree, provided they have 
roots in a circle and all other roots are far away from this circle. 

It is useful to have a decision criterion for the existence of factors of 
degree greater than 1. We will do that for degree-2 factors, since this is the 
interesting case for real polynomials. See Ostrowskii |33[ for more results. 

Clearly, it is enough to consider the case of a polynomial / of degree 2: 



f(x) = f2X 2 + flX + f = X 2 + (-& " + 0C2 

In case |0| 10 1 5 we have: 



log 



I./2 



+ log 



l./i 



In case R = |(i 



\fi\J Vl/ol 
|0|, we can bound: 



>log|^»0 

I (.2 1 



(28) 



(29) 



log (H) +iog ([f) <log4 



(30) 



However, we should look at the renormalized Newton diagram. In that 
case, we should divide equations (P§Q and (^) by 2 k , and expect that if 
10 1 7^ 101 then the following holds: 



> a 



(31) 
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Zeros of a generic polynomial of degree 1000 




-25 -20 -15 -10 -5 



10 15 



Figure 2: Zeros of a random polynomial 



where a is an a priori bound on root moduli separation. It may be obtained 
from a probabilistic analysis (Lemma |6|), or from any other a priori knowledge 
on the polynomial; for the choice of a, notice that if |Ci| = \C2\1 



-" '" K, 7T7 

.l/o| 



< 2- fc log4 



(32) 



Thus we may choose k such that 2~ fe log4 < a. In case equation (|3lD is 
not satisfied, it is reasonable to assume that the two roots have the same 
modulus indeed. 



6 Numerical Results 

The results discussed below are tentative, and our algorithm deserves further 
experimentation. Moreover, at this moment we are using a tentative algo- 
rithm to find the argument of the roots. Those matters will be dealt with 
in a subsequent paper (PSfl). The purpose of this section is to illustrate 
the numerical properties of Renormalized Graeffe Iteration, when applied to 
random real polynomials. 

A C implementation of our algorithm was tested for pseudo-random real 
and complex polynomials. The arguments of the solutions were recovered 
using an algorithm derived from the Renormalized Graeffe Iteration. 
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IMPLEMENTATION DETAILS: Polynomials with zero coefficients do 
arise in practice, and we explain now how to deal with the 'log of zero' 
problem. The IEEE floating point arithmetic, implemented by most modern 
computers, has a few useful features for this sort of situation. When one 
asks a computer to produce logO, it returns a special IEEE value called — oo. 
This is not an error, and further calculations can be carried out (as long as 
they are defined). If they are not defined (e.g. oo — oo) then another special 
value, called a 'not a number' is returned. 

This last case (oo — oo) can appear during the computation of a renor- 
malized sum. It can be dealt by testing a and b in Example |^ for finitesness. 
In case a or b is infinite, then c should be set to oo. For an introduction of 
IEEE arithmetic, see jBJ pages 45-48 or \T2\ pages 9-15. The correctness 



of the results was certified by estimates as in p6 |. 

Experiments were performed in a Pentium 66 system running Linux op- 
erating system. Since the objective here was to illustrate the asymptotic be- 
haviour of the algorithm, we did not perform experiments in other systems. 
Those would be necessary if one wanted to compare with other algorithms 
with same asymptotic properties. However, this goes far beyond the scope 
of this paper. 

The table below shows the average and median user time in a Pentium- 
based computer, using 'long double' precision. Time does not include val- 
idation time. Ten pseudo-random polynomials were tested for each degree. 
The actual experimental data is plotted in Figure 0. 



Degree 


Real polynomials 


Complex polynomials 




avg time (s) 


median time(s) 


avg time (s) 


median time (s) 


100 


0.87 


0.87 


1.19 


1.18 


200 


3.23 


3.24 


4.35 


4.30 


300 


7.07 


7.07 


9.99 


9.73 


400 


12.60 


12.58 


17.28 


17.10 


500 


19.41 


19.38 


26.75 


26.71 


600 


27.67 


27.73 


37.02 


35.96 


700 


37.50 


37.32 


51.30 


49.91 


800 


48.89 


48.72 


65.39 


63.61 


900 


61.56 


61.06 


79.28 


78.74 


1000 


75.89 


75.47 


102.33 


101.80 



We also computed (approximately) the relative separation of the moduli 



32 



Time (s) for solving real polynomials 




1000 



Figure 3: Timing for 100 random real polynomials 



of the solutions Q: 

\Ci\ - 10 



mm 



12 



KiKioi VTT]Q\ 

The values obtained are also plotted in figures | to |. 

Further experimentation is necessary to obtain data about polynomials 
of degree ^> 1000. Indeed, due to underflow, we cannot represent random 
high-degree polynomials in the usual floating point representation. 
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Time (s) for solving complex polynomials 



100 - 




1000 



Figure 4: Timing for 100 random complex polynomials 
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Moduli separation of real polynomials 



0.001 - 



0.0001 



1e-05 - 



1e-06 



1e-07 



1000 



Figure 5: Separation for 100 random real polynomials 
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Moduli separation of complex polynomials 



0.0001 - 



1e-06 - 




1000 



Figure 6: Separation for 100 random complex polynomials 
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